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Summary 

The analysis of point patterns with linear struetures is of interest in many applieations. To de- 
teet anisotropy in sueh eases, in partieular in ease of a eolumnar strueture, we introduee a fune- 
tional summary statistie, the eylindrieal JT-funetion, whieh is a direetional iT-funetion whose 
strueturing element is a eylinder. Further we introduee a elass of anisotropie Cox point pro- 
eesses, ealled Poisson line eluster point proeesses. The points of sueh a proeess are random 
displaeements of Poisson point proeesses defined on the lines of a Poisson line proeess. Parame¬ 
ter estimation based on moment methods or Bayesian inferenee for this model is diseussed when 
the underlying Poisson line proeess is latent. To illustrate the methodologies, we analyze two- 
and three-dimensional point pattern data sets. The three-dimensional data set is of partieular in¬ 
terest as it relates to the minieolumn hypothesis in neuroseienee, elaiming that pyramidal and 
other brain eells have a eolumnar arrangement perpendieular to the surfaee of the brain. 

Some key words: Anisotropy; Bayesian inference; Directional if-function; Minicolumn hypothesis; Poisson line pro¬ 
cess; Three-dimensional point pattern analysis. 


1. Introduction 

Frequently in the spatial point proeess literature, isotropy, i.e., distributional invarianee under 
rotations about a fixed loeafion in spaee, is assumed, fhough if is offen unrealisfie. Anisofropy 
of spatial poinf proeesses has usually been sfudied by summarizing fhe information of observed 
pairs of poinfs, ineluding fhe use of direefional iC-funefions or relafed densifies (Ohser & Sfoyan, 
1981; Sfoyan & Benes, 1991; Sfoyan, 1991; Sfoyan & Sfoyan, 1995; Guan ef ah, 2006; Illian 
ef ah, 2008; Redenbaeh ef ah, 2009), speefral and wavelef mefhods (Mugglesfone & Renshaw, 
1996; Rosenberg, 2004; Nieolis, Mafeu & D’Ereole, 2010), and geomefrie anisofropie pair eor- 
relafion funefions (Mpller & Toffaker, 2014). The applieations eonsidered in fhese referenees 
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Figure 1. Poisson line cluster process: Left panel: A sim¬ 
ulated realization of a Poisson line cluster point process 
within a three-dimensional box. The realizations of the 
Poisson line process (solid lines) and the Poisson point 
processes on the lines (circles) are also shown. The dotted 
lines indicate how the points on the lines have been dis¬ 
placed to new positions (filled circles) and they specify the 
clusters. Right panel: The same simulated realization of a 
Poisson line cluster point process and different choices of 
cylinders centered at different points of the process. 


except Redenbach et al. (2009) and Illian et al. (2008) are for two-dimensional but not three- 
dimensional point patterns. 

There are point patterns where points lie approximately along straight lines, cf. the examples 
of applications and references in Mpller & Waagepetersen (2016), called point patterns with 
linear structures. This paper focuses on detecting and modelling such point patterns observed 
within a bounded subset of IR'^, d> 2, where the cases d = 2 and d = 3 are of main interest. In 
particular we study columnar structures. 

Section 2 introduces the cylindrical iT-function, a directional iT-function whose structuring 
element is a cylinder which is suitable for detecting anisotropy caused by columnar or other 
linear structures in spatial point patterns. This is an adapted version of the space-time iT-function 
(Diggle et al., 1995; Gabriel & Diggle, 2009). Section 3 concerns a new class of point processes, 
Poisson line cluster point processes, whose points cluster around a Poisson line process. The 
left panel in Figure 1 illustrates how such a process is constructed: lines are generated from an 
anisotropic Poisson line process, independent stationary Poisson point processes are generated 
on the lines, and their points are randomly displaced, resulting in the Poisson line cluster point 
process. We consider the Poisson lines and the points on the lines as latent, so the clusters of 
the Poisson line cluster point process are also hidden. Section 3 also discusses a moment based 
approach and a simulation-based Bayesian approach for inference, where in the latter case we 
estimate both the parameters of the model and the missing lines. 

Sections 2-3 apply our methodology to the data sets in Figure 2. The left panel shows a two- 
dimensional point pattern data set recorded by Mugglestone & Renshaw (1996), namely the 
locations of 110 chapels in the Welsh Valleys, United Kingdom, where the clear linear orientation 
is caused by four more or less parallel valleys. For a three-dimensional point pattern data set it 
is often difficult to detect anisotropy by eye, but the cylindrical iT-function will be useful, as 
illustrated later in connection to the right panel, which shows the locations of 623 pyramidal cells 
from the Brodmann area 4 of the grey matter of the human brain collected by the neuroscientists 
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Figure 2. Data sets: Left panel: Locations of 110 chapels 
in Wales, United Kingdom, observed in a square win¬ 
dow (normalized to a unit square). Right panel: Nucleolus 
of 623 pyramidal cells in an observation window of size 
508 X 138 X 320 firrP. 


at the Center for Stoehastie Geometry and Bioimaging, Denmark. Aeeording to the minieolumn 
hypothesis (Mounteastle, 1957), brain eells, mainly pyramidal eells, should have a columnar 
arrangement perpendicular to the pial surface of the brain, i.e., a columnar arrangement parallel 
to the x^^^-axis indicated in Figure 2, and this should be highly pronounced in Brodmann area 
4. However, this hypothesis has been much debated, see Rafati et al. (2016) and the references 
therein. 

We use the chapel data set mainly for illustrative purposes and for comparison with previous 
work. Investigations of the minieolumn hypothesis have so far only been done in two dimensions 
except for the three-dimensional analysis in Rafati et al. (2016). The present paper details the 
methodology and provides a more thorough analysis of the pyramidal cell data set, shedding 
further light on the validity of the minieolumn hypothesis. 

Throughout Sections 2-3 we assume stationarity. Section 4 discusses the choice of a cylinder 
as the structuring element of the iT-function and extensions of this function and the Poisson line 
cluster point process in a non-stationary setting. 


2. The cylindrical TT-function 
2-1. Setting 

Throughout this paper we make the following assumptions and use the following notation. 

Unless otherwise stated, we consider a stationary point process X defined on IR'^, with finite 
and positive intensity p, and where we view 26 as a locally finite random subset of Here 
stationarity means that the distribution of X is invariant under translations in R'^, and p\B\ is 
the mean number of points from X falling in any Borel set H C R'^ of volume \B\. We assume 
that X has a pair correlation function g{x) defined for all x € R*^. Intuitively, if xi,X 2 G R'^ are 
distinct locations and Bi,B 2 are infinitesimally small sets of volumes dxi, dx 2 and containing 
xi,X 2 , respectively, then p^g{xi — X 2 ) dxi dx 2 is the probability for X having a point in each 
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of Bi and B 2 . For further details on spatial point proeesses, see Mpller & Waagepetersen (2004) 
and the referenees therein. 

We view any veetor X = ..., G IR*^ as a eolumn veetor, and ||x|| = + 

... + as its length. For ease of presentation, we assume that no pair of distinet points 

{xi, X 2 } C X is sueh that u = (xi — X 2 )/||xi — X 2 II is perpendieular to the x'^'^^-axis. This will 
happen with probability one for the models eonsidered later in this paper. 

Let = {u = {u^^\... G : ||tt|| = 1} be the unit sphere in R*^ and = 

(0,..., 0,1) its top point. Denote o = (0,..., 0) the origin of R*^. Consider the d-dimensional 
cylinder with midpoint o, radius r > 0, height 2t > 0, and direction e^: 

C{r,t) = |x = (x(i),...,x('^)) G R*^; + • • • + < r^, < f} . 

For u G denoting Ou an arbitrary d x d rotation matrix such that u = Oued, then 

Cuir,t) = OuC{r,t) 

is the d-dimensional cylinder with midpoint o, radius r, height 2t, and direction u. 

2-2. The cylindrical K-function 

Recall that the second order reduced moment measure K, with structuring element B C R*^, a 
bounded Borel set, is given by 

IC{B) = j g{x) dx 
Jb 

(Mpller & Waagepetersen, 2004, Section 4.1.2). Ripley’s iT-function (Ripley, 1976, 1977) is 
obtained when 77 is a ball and it is not informative about any kind of anisotropy in a spatial point 
pattern. 

To detect preferred directions of linear structures in a spatial point pattern, in particular a 
columnar structure, we propose a cylinder as the structuring element and define the cylindrical 
TT-function in the direction u by 

Ku{r,t) = / g{x)dx, u£S^~^, r>0,t>0. (1) 

J C„(r,t) 

Intuitively, pKu {r, t) is the mean number of further points in X within the cylinder with midpoint 
at the typical point of X, radius r, and height 2t in the direction u. For example, a stationary 
Poisson process is isotropic, has g = 1, and 

Ku{r,t) = 

where oJd-i = -|- l)/2} is the volume of the (d — 1)-dimensional unit ball. For 

(7 = 3, 7T(o,o,i) is similar to the space-time TT-function in Diggle et al. (1995) and Gabriel & 
Diggle (2009), when considering the x^^i-axis as time and the (x^^i, xi^i)-plane as space. 

If FF C R*^ is an arbitrary Borel set with 0 < |FF| < 00 , then by standard methods (Mpller & 
Waagepetersen, 2004, Section 4.1.2) 

Ku{r,t)= 21 WI l{xi G FF,X 2 - xi G C'„(r,f)} , r > 0, f > 0, (2) 

X\^X2^X\X\^X2 

where 1 denotes the indicator function, and by stationarity the right-hand side does not depend 
on the choice of W. This provides a more general definition of K^, since (2) does not require 
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the existence of the pair correlation function. Equation (2) becomes useful when deriving non- 
parametric estimates in Section 2-3. 


2-3. Non-parametric estimation 

Given a bounded observation window IE C IR*^ and an observed point pattern {rci,..., Xn} C 
W with n > 2 points, we consider non-parametric estimates of the form 


Ku{r, t) 


i ^ vJu{xi, Xj)l{xj - Xi e Cu{r, t)}. 


(3) 


Here is a non-parametric estimate of and Wu is an edge correction factor. If X is isotropic, 
Ku{r, t) does not depend on u and this should affect the choice of Ku{r, t). On the other hand, 
as illustrated in the right panel of Figure 1 and in Section 24, to detect a preferred direction 
of linearity in a spatial point pattern, we suggest using an elongated cylinder with t > r and 
considering different directions u. Then we expect a largest value of Ku{r,t) to indicate the 
preferred direction, but a careful choice of r and t may be crucial, cf. Section 4. Furthermore, 
since Ku = K-u, we need only to consider the case where u = {u^^\ ..., is on the upper 
unit-sphere, i.e., uW > 0 . 

Specifically, we use = ^{n — l)/|VFp, see e.g. Illian et al. (2008), and the translation cor¬ 
rection factor (Ohser & Stoyan, 1981) 


Wu{xi,X 2 ) = l/\w (4) 

where Wx denotes translation of the set IE by a vector x G IR*^. Then, by Femma 4.2 in Mpller 
& Waagepetersen (2004), if p^ is replaced by p^ in (3), we have an unbiased estimate of Ku- As 
in Figures 1-2, if IE is rectangular with sides parallel to the axes and of lengths ai,..., > 0, 

d 

|1E n Wx,-x,l = xi,X2 G w, 

i=l 


where denotes the f’th coordinate of Xj {j = 1 , 2 ). 

For d = 3, IE = [0, oi] x [0, 02 ] x [0, 03 ], and u = (0,0,1), another choice is a combined 
correction factor 


W(0,0,1){X1,X2) = 


1 + 1 2x 


M (3) 


— Xo 


0 [0,03]) 


03 at — 


41 ) ^( 1 ) 


- X 


02 - 


„( 2 ) ..( 2 ) 
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where the numerator is a temporal correction factor and the denominator is the reciprocal of a 
spatial correction factor; similarly we construct combined correction factors when u = ( 1 , 0 , 0 ) 
or n = (0,1,0). Instead of this spatial correction factor, which is of a form similar to (4), Diggle 
et al. (1995) used an isotropic correction factor, but this is only appropriate if X is isotropic in 
the (xi, X 2 )-plane. The temporal correction factor is the same as that used in Diggle et al. (1995). 

We prefer the translation correction factor (4), since this does not restrict the shape of IE and 
the choice of u. In a simulation study with d = 3, IE = [0,1]^, and X a Poisson line cluster point 
process as defined in Secfion 3, we obfained similar resulfs when using fhe franslafion and fhe 
combined correcfion faclors. 
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Figure 3. Chapel data set: The left panel shows the 
non-parametric estimate -^{cos ip.sirK^a) (f, f) versus ip for 
four different combinations of r and t with the curves 
from the top to the bottom corresponding to (r, t) — 
(0.2,0.4), (0.2,0.3), (0.1, 0.3), (0.1, 0.2) and the two 
vertical lines corresponding to 113° and 124°. The middle 
panel shows Tf{cos<, 5 ,sin(, 3 )(f, t) versus r for different val¬ 
ues of ip and with t — 0.3 with the solid curves from the 
top to the bottom corresponding to 20°, 45°, 170°, and the 
dotted curve corresponding to = 117°. The right panel 
shows a non-parametric estimate of the point pair orienta¬ 
tion distrihution function with ri = 0.05 and r 2 = 0.15. 

For more details, see Section 2-4. 


2-4. Examples 

Non-parametric estimates of the cylindrical iF-function for the two-dimensional chapel data 
set and the three-dimensional pyramidal cell data set are shown in Figures 3 and 4, respectively. 
Below we comment on these plots. Further examples are given in Rafati et al. (2016). 

To detect the main direction in the chapel point pattern, the left panel in Figure 3 shows, 
for four different combinations of r and t, i.e., (0.1,0.2), (0.1, 0.3), (0.2, 0.3) and (0.2,0.4), 
plots of Ku{r,t) versus p, where rt = (cos((/9), sin((^)). These four curves are approxi¬ 
mately parallel, and a similar behaviour for other choices of r = 0.05,0.1,0.15,0.2 and t = 
0.15,0.2,0.25,0.3,0.35,0.4 with t > 2r was observed. In a previous analysis, Mpller & Tof- 
taker (2014) estimated the orientation of the chapel point pattern to be between 113° and 124°. 
This interval, which is specified by the vertical lines in the left panel of Figure 3, is in close 
agreement with the maximum of the Ku{r, f)-curves. The middle panel in Figure 3 shows plots 
of Ku{r, t) versus r when t = 0.3. For the dotted curve in the middle panel, tp = 117° is the av¬ 
erage of the four maximum points of ip corresponding to the four curves in the left panel, while 
for the three other curves in the middle panel, values of ip not included in the interval [113°, 124°] 
have been chosen. The clear difference between the dotted curve and the other curves indicates 
a preferred direction in the point pattern which is about 117°. This is also confirmed by fhe 
righf panel in Figure 3, which shows a non-paramefric esfimafe of fhe poinf pair orienfafion 
disfribufion funclion given by equafion (14.53) in Sfoyan & Sfoyan (1995) and implemented in 
spat St at (Baddeley & Turner, 2005). This is a kernel esfimafe which considers fhe direcfion 
for each pair of observed poinfs fhaf lie more fhan ri = 0.05 and less fhan r 2 = 0.15 unite aparf. 
For ofher values of r 2 < 0.27 we reached similar conclusions, buf for higher values of r 2 fhe pair 
orienfafion disfribufion function did nof show a clear preferred direcfion in fhe dafa. 
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Figure 4. Pyramidal cell data set: Non-parametric estimates 
Ku{r, 80) — IGOyrr^ versus r when t = 80 and the cylin¬ 
der is along the -axis (dotted line), the -axis (dotted 
line), or the x^^^-axis (solid line). The grey region speci¬ 
fies a 95% simultaneous rank envelope computed from 999 
simulations under complete spatial randomness. For more 
details, see Section 2-4. 


By the minicolumn hypothesis, the pyramidal cell data set has a columnar arrangement 
in the direction of the x^^^-axis indicated in Figure 2 (Rafati et ah, 2016). Figure 4 shows 
that the cylindrical iF-function is able to detect this kind of anisotropy: The three curves are 
Ku{r, 80) — IGOvrr^ for 0 < r < 20 and u parallel to one of the three main axes, where lOOvrr^ 
is the value of Ku{r, 80) under complete spatial randomness, i.e. a stationary Poisson point pro¬ 
cess model; we only make this comparison in order to see any deviations from complete spatial 
randomness. The solid curve corresponding ton = (0, 0,1), i.e., when the direction of the cylin¬ 
der is along the x^^^^-axis, is clearly different from the two other cases where tt = (1, 0, 0) or 
ii = (0,1,0). The grey region is a so-called 95% simultaneous rank envelope (Myllymaki et ah, 
2016) obtained from 999 simulated realizations under complete spatial randomness; Myllymaki 
et al. (2016) recommended 2499 simulations, however, for the cylindrical iT-function considered 
in Figure 4, 999 simulations seemed sufficient since results were produced that were similar to 
those using 2499 simulations. Roughly speaking, under complete spatial randomness, each of the 
estimated cylindrical iT-functions is expected to be within the grey region with estimated prob¬ 
ability 95%. While the curves for u = (1,0,0) and u = (0,1, 0) are completely within the grey 
region, the curve for u = (0,0,1) is clearly outside for a large range of r-values. In fact, for the 
null hypothesis of complete spatial randomness, considering the rank envelope test (Myllymaki 
et ah, 2016) based on Ku{r, 80) when 0 < r < 20 and u = (0,0,1), the p-value is estimated to 
be between 0.1% and 0.18%, showing a clear deviation from the null hypothesis of complete 
spatial randomness. 

These examples illustrate that the cylindrical iT-function is a useful functional summary statis¬ 
tic for detecting preferred directions and columnar structures in a spatial point pattern. Particu¬ 
larly, for the pyramidal cell data set and in accordance to the minicolumn hypothesis, there is a 
pronounced columnar arrangement in the direction of the x^^^-axis. 
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3. The Poisson line cluster point process 
3T. Definition of Poisson line cluster point processes 

Motivated by the analysis in Seetion 2-4, in partieular that of the pyramidal eell data set, 
we now introduee a model for a point proeess X with eolumnar strueture. It is a hierarehieal 
eonstruetion, using various latent proeesses speeified briefly in (Al)-(Cl) and later in more detail 
in (A2)-(C2), while X is given in (Dl) and (D2) below; the left panel of Figure 1 is helpful in 
this eonneetion: It eonsists of generating: 

(Al) a Poisson line proeess L = • • •} of lines li, i.e., infinite, direeted, and straight lines; 

only parts of these lines are shown in Figure 1; 

(Bl) on eaeh line li, a Poisson proeess Qi illustrated by unfilled points in Figure 1; 

(Cl) a new point proeess Xi obtained by random displaeements in IR'^ of the points in Qi illus¬ 
trated by filled points in Figure 1; 

(Dl) finally, X as the superposition of all the X*. 

Then we eall X a Poisson line eluster point proeess, sinee its points eluster around the Poisson 
lines, and we eall eaeh Xi a eluster. 

In eonneetion to the more detailed eonditions (A2)-(D2) we need the following notation. 
Let • denote the usual inner produet on IR*^. For u G let ru*- = {x G R'^ : x • u = 0} be the 
hyperplane perpendieular to u and eontaining o, A„x the {d — 1)-dimensional Lebesgue measure 
on n-*-, and Pu^{x) = x — (x • u)u the orthogonal projeetion of x G R'^ onto Let H = ejp 
and A = A x, i.e., H is the hyperplane perpendieular to the x^'^^-axis. Further, let /c be a density 

funetion with respeet to Lebesgue measure on As in Seetion 2-1, suppose we have speeified 
for eaeh uGS'^^^adxd rotation matrix Ou sueh that u = OuGd- We then define a density with 
respeet to A^x by 

k^±{Ou{x^^\ .. ■ = k{x^^\ ... (x«,...,x(‘^-i)) G R-^-^ 

In other words, when eonsidering eoordinates with respeet to the d — 1 first eolumns in Ou, the 
distribution under ku± is the same as under k. Furthermore, for the line proeess L we use the 
so-ealled phase representation (Chiu et al., 2013) and assume that with probability one, L has 
no line eontained in H. Thereby a line I = l{y,u) in L eorresponds to its direetion u G 
and its interseetion point y in H. Thus L = {h,l 2 , ■ ■ •} can be identified by a point proeess 
4> = {(yi,ui), ( 2 / 2 , 112 ), ...} C H X sueh that k = l{yi,Ui) (i = 1, 2,...) and 4> C iT x 
(S<t-i almost surely. 

In addition to (Al)-(Dl) we assume that: 

(A2) is a Poisson proeess with intensity measure l3X{dy)M{du), where /3 > 0 is a parameter 
and M is a probability measure on deseribing the direetion of a typieal line, where 
^(S^-i Pi — Q. assume that f M{du) < 00 , where is the last eoor- 

dinate of u; this assumption will be needed when we later in (7) speeify the intensity and 
rose of direetions of the line proeess; 

(B2) eonditional on <1>, we have that Qi,Q 2 , ■ ■ ■ are independent stationary Poisson proeesses 
on (i, ( 2 , • • •, respeetively, with the same intensity a > 0; 

(C2) eonditional on <I> and Qi, Q 2 , • • •, we have that Xi, X 2 ,... are independent point proeesses 
and eaeh Xj is obtained by independent and identieally distributed random displaeements 
of the points in Qi following the density /c^x; thus X* is a Poisson proeess on R^ with 
intensity funetion 

Aj(x) = a/c„x{p„x(x - 2 /j)}, X G R"*; 

i i 


(5) 
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(D2) hence the superposition X = is a Cox process driven hy A = ^ ■ Aj, i.e. X con¬ 

ditional on <I> is a Poisson process with intensity function 

OO 

A{x) = a^k^±{p^±{x-Vi)}, X G R"*. (6) 

i=\ 

Some comments are in order. 

The processes L, A, and X are stationary, the distribution of L is given hy (/3, M), and the 
distribution of X is determined by (/3, M, a, k). 

By (C2), conditional on L, for each line k ^ L and each point qij G Qi, there is a corre¬ 
sponding point Xij G Xi such that the random shift Zij = Xij — qij follows the density . We 
could have defined the Poisson line cluster point process by letting the displacements follow a 
distribution on IR'^ rather than a hyperplane, or more precisely by letting Zij follow a density 

— Puiizij)} 1 Zij G IR , 

where is a density function with respect to Lebesgue measure on the line k — yt = {tui : 
f G IR}. However, since the part of the displacements running along the line k just corresponds to 
independent displacements of a stationary Poisson process, this will just result in a new stationary 
Poisson process with the same intensity, see, e.g.. Section 3.3.1 in Mpller & Waagepetersen 
(2004), and so there is essentially no difference. 

The fact in (D2) that is a Cox process becomes important for the calculations and the 
statistical methodology considered later in this paper. 


3-2. Intensity and rose of directions for the Poisson line process 
We have specified fhe disfribufion of fhe Poisson line process L by (/?, M). This is useful for 
compufafional reasons, buf when interpreting resulfs if is usually more nafural fo consider fhe 
infensify and fhe rose of direcfions of L, which we denofe hy pi and TZ, respectively. Formal 
definitions of fhese concepfs are given in Appendix A, where if is shown fhaf for any Borel sef 
B C 


PL = p j l/\u^<^'l\M{du), n{B) = j l/\u^‘^'>\M{du)^ y'l/|u("')|M(du). (7) 

In words, pi is the mean length of lines in L within any region of unit volume in IR'^, and TZ is 
the distribution of the direction of a typical line in L, see, e.g., Chiu et al. (2013). 

Equation (7) establishes a one-to-one correspondence between {pL,TZ) and (/3, M), where 


/3 = PL 


u^^^\IZ{du), 



u^^^\n{du). 


( 8 ) 


Consequently, we can choose pi as any positive and finite parameter, and TZ as any probability 
measure on Moreover, (3 < pl where the equality only holds when TZ is concentrated with 
probability one at ie^. We call this special case the degenerate Poisson line cluster point process. 

For the rose of directions, we later use a von Mises-Fisher distribution with concentration 
parameter k > 0 and mean direction p, G This has a density /(• | p,k) with respect to the 
surface measure on 

/(ii I = Crf(K)exp(Kp-n), Cd{n) = . -—, u G S"* 


(9) 
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where denotes the modified Bessel funetion of the first kind and order d. Note that L and X 
are then isotropie if and only if k = 0, in whieh ease the ehoiee of ^ plays no role. For k > 0, 
the direetions of the lines in L are eoneentrated around fi, and so the elusters in X have preferred 
direction /x. When /x = zhe^, in the limit as k — )■ oo, we obtain the degenerate Poisson line cluster 
point process. 


3-3. Finite versions of the Poisson line cluster point process and simulation 
Suppose we want to simulate the Poisson line cluster point process within a bounded region 
W C IR'^, i.e. the restriction Xw = ^ n W. Then we need a finite approximation of which 
will also be used when we later discuss Bayesian inference, as follows. Consider a bounded 
region IFext 2 W and let 

S = {{y,u) e H X : l{y,u) n Wext / 0} 

be the set of all lines hitting IFext- We want to choose IFext as small as possible but so that it is 
very unlikely that for some line U £ L with {yi,Ui) 0 S, Xi has a point in W. Then our finite 
approximation is ‘hs = n 5, and (i) we simulate <I> 5 , and (ii) conditional on <I >5 we make an 
approximate simulation of Xy/ as a Poisson process with intensity function 

kw{x) = a ^ k^v{p^v{x - y)}, x G W, (10) 


cf. ( 6 ). We detail (i)-(ii) below. 

Here (ii) is rather straightforward: Suppose we have simulated ‘hs' and consider any {yi, Ui) G 
The projection of W onto k is the bounded set lw,i = {x G U ■. {x + uf) n VF / 0}. In 
accordance to (B2), we simulate a Poisson process TV,* with intensity a on ly/,i- Displacing 
the points in Yjy j as described in (C2) we obtain a Poisson process Xy/^i with intensity func¬ 
tion (5) but restricted to \Jx&i^.{x + xx-*-). The approximate simulation of Xyy is then given by 

In (i) we assume for simplicity and specificity that TZ follows the von Mises-Fisher density 
(9). Denote the surface measure on Then ( 8 ) implies that 

/3A(dy)/x(d;u) = pifiy \ p, K)\u^'^'^\X{dy)ud-i{du), 

i.e., is a Poisson process on S with intensity function 

x{y,u I pl,P,k) = I p,k) 

with respect to the measure A(dy)r'(X_i(du). First, we therefore simulate the Poisson distributed 
counts with mean pLl{p, n) where 

/(/x,k) = j \u^'^^\f{u I p, K)dX{dy)ud-iidu) = j A(J„)/(xx | p, K)i'd-i{du) 

where Ju = {y G FI : l{y, u) n IFext / 0}- Second, we simulate each (y, u) G <hs’ with density 
proportional to |xx*^'^^|/(xx | p, n) for y G Ju and zero otherwise. Here we use rejection sampling. 

For example, if d = 2 and IFext = [—ci,a]‘^ is a square centered at the origin, then for u = 
(cos p, sin p), we have Ju = Jip x { 0 } with 

J _ f [—a cot p — a, a cot p + a] , 0 < p < 7 r /2 or tt < p < ‘dtt 
^ I [a cot (/? — a, a — a cot , 7 r /2 < (^ < vr or 37 r /2 < < 27r. 


(11) 
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Further, 


A(Jn) 


2a + 2a cot p, 0 < (p < 7r/2 or tt < (p < 37r/2, 
2a — 2a cot p, 7r/2 < p < tt or 37r/2 < p < 27r, 


( 12 ) 


and 


/* 7 r/Z PIT 

I{p,,K)=2a / (sin(^ + cos(/?)/(u I/r, k) d(/? + 2a / (sin(/9 — cos(/9)/(rt |//, k) d(/9 

JO A/2 

r 3 nj 2 p 27 t 

— 2a {sin pcos p)f{u \ p, k) dp — 2a / (sine/? — cos 99)/(rt | ^, k) d(^, 

J TT J37T/2 


(13) 


which can be evaluated by numerical methods. Furthermore, for p = (cos 6*, sin 0) and y = 
(y(i)^y(2)), unnormalized density \u I p,k) = G J(p)| sin(/9| exp{Kcos((^— 

0)} is just with respect to Lebesgue measure dy^^^ d(/? on IR x [0, 27r). Finally, when doing rejec¬ 
tion sampling, we propose p from /(• | p, k) and y^^^ from the uniform distribution on J^p, and 
accept {p, y^^)) with probability | sin(/9|. 


3-4. Moments of the Poisson line cluster point process 
Since X is a Cox process with driving random intensity A, moment properties of X are deter¬ 
mined by moment properties of A. This section focuses on first and second order moments. 

The Poisson line cluster point process X has intensity p and pair correlation function g 

P = ^{A-(o)}, p‘^g{x) = F;{A(o)A(x)}, X G R"*. (14) 

Appendix B verifies fhaf 

P = apL (15) 

and 

g{x) = 1 + — ( * fc„±{p„±(x)}7^(du), X G R"*, (16) 

PL J 

where k^^{p^^{x)} = k^^{—p^^{x)} and * denotes convolution, i.e., 

kuX * ~kux{Pu^{x)] = j k^r{p^±{x) - y}Kx{y) A„±(dy). 

Thus y > 1, reflecting fhe clusfering of fhe Poisson line clusfer poinf process. Evaluation of fhe 
integral in (16) may require numerical mefhods. For example, if k{-) = /(• | a^) is fhe densify 
of fhe {d — 1)-dimensional zero-mean isofropic normal disfribufion wifh variance cr^ > 0, fhen 

V *kux{Pux{x)} =exp{-||p„x(x)||V ( 4 ^ 72 )}/( 47 r^T 2 )^‘^■^^/^ (17) 

3-5. Moment based inference 

The likelihood for a paramefric Poisson line cluster poinf process model is complicafed be¬ 
cause of fhe hidden line process and fhe hidden poinf processes on fhe lines, fhough if can be 
approximated using a missing dafa Markov chain Monfe Carlo approach, see, e.g., Mpller & 
Waagepefersen (2004). A Bayesian Markov chain Monfe Carlo approach is used in Seefion 3-6 
where fhe missing dafa is included info fhe posferior. Simpler procedures for paramefer esfi- 
mafion are composite likelihood (Guan, 2006; Mpller & Waagepefersen, 2007) and minimum 
confrasf mefhods (Diggle & Grafton, 1984) based on (15)-(16). Since g is hard to compute in 
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general, this seetion foeuses on sueh simple proeedures in the speeial ease of a degenerate Pois¬ 
son line eluster point proeess whieh e.g. is a relevant model for the pyramidal eell data set shown 
in Figure 2. For speeifieity we assume as in (17) that k{-) = f{- \ Then the unknown param¬ 
eters are /3 = PL > 0, a > 0, and > 0. 

Suppose a realization of Xw is observed within a region of the produet form W = D x I, 
where D C and / C IR are bounded sets. To estimate the unknown parameters we notiee 
the following. Let Xj denote the projeetion of X n x I) onto H. Sinee we eonsider a 

degenerate Poisson line eluster point proeess, the x^'^^-eoordinates of the points in Xw are in¬ 
dependent and identieally distributed uniform points on I whieh are independent of Xj. Thus 
X] is a suffieient statistie for a, cj^). Note that Xj is a Cox proeess driven by the random 
intensity funetion 

OO 

T{x) = a\I\'^f{x-yi\a‘^), x e H, 

i=l 

where <1> = {(yi, e^), (y2, erf); • • •} can be identified by the stationary Poisson proeess 
{ 2 / 1 ) 2/21 •• •} on H with intensity pi- Therefore X/ is a modified Thomas proeess (Mpller & 
Waagepefersen, 2004) wifh infensify pj = p\I\ and by (16)-(17) pair eorrelafion funefion 

Paramefer esfimafion based on [pj, gj) and using a eomposife likelihood or a minimum eon- 
frasf mefhod is slraighfforward (Mpller & Waagepefersen, 2007). Then, when eheeking a fiffed 
Thomas proeess, we should nol reuse fhe infensify and fhe pair eorrelafion funefion. Below we 
use instead fhe funelional summary sfafisfies fhe empfy spaee funefion F, fhe nearesf-neighbour 
funefion G, and fhe J-funefion (Mpller & Waagepefersen, 2004). 

As an example, for fhe fhree-dimensional pyramidal eell dafa sef in Figure 2 in aeeordanee 
wifh fhe minieolumn hypofhesis, we eonsider a degenerate Poisson line elusfer poinf proeess. 
This has a eolumnar arrangemenf in fhe direefion of fhe x^^^-axis and fhe observafion window 
is of fhe same form as deseribed above wifh D = [0, 508] x [0,138] and / = [0,320]. The firsl 
panel in Figure 5 shows fhe empirieal eumulafive disfribufion funefion of fhe x^^^-eoordinafes of 
fhe pyramidal eell poinf paffern dafa sef; fhere is no elear indieafion of a deviation from a uniform 
disfribufion, in agreemenf wifh our sfafionarify assumpfion. 

When tiffing fhe modified Thomas proeess for fhe projeefed poinf paffern onto D, for bofh 
eomposife likelihood and minimum eonfrasf estimation, we used fhe spatstat (Baddeley 
& Turner, 2005) funefion kppm. We obfained fhe minimum eonfrasf esfimafes pi = 0.024, 
a = 0.37/320 = 0.0012, and = 15.04, and similar estimates were obfained by a eomposife 
likelihood mefhod. The fhree lasf panels in Figure 5 show fhe non-paramefriely esfimafed F, G, 
and J-funelions for fhe projeefed pyramidal eell poinf paffern onfo D, fogefher wifh 95% simul- 
faneous rank envelopes eompufed from 4999 simulated poinf pafferns under fhe fiffed Thomas 
process. The p-values for fhe rank envelope fesf (Myllymaki el ah, 2016) for G, F, and J- 
funclions are wifhin fhe intervals [0.851, 0.852], [0.732, 0.733], and [0.623,0.625], respectively, 
providing no evidence againsf fhe filled model. 

3-6. Bayesian inference 

Suppose a non-emply realization Xy/ = {xi,... ,Xn} is our dafa, where 4F C R'^ is a 
bounded observafion window, and we model X as a Poisson line elusfer poinf process wifh 
k^± (y) = /(y I cr^) and TZ following fhe von Mises-Fisher density f{u \ p, n) given by (9). As 
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Figure 5. Summary statistics for the pyramidal cell point 
pattern data set: Empirical cumulative distribution function 
of the a:-coordinates of the pyramidal cell point pattern 
data set (top left), and non-parametric estimates of G (top 
right), F (bottom left), and J (bottom right) for the pro¬ 
jected pyramidal cell point pattern onto D (solid lines), 
together with 95% simultaneous rank envelopes (gray re¬ 
gions) computed from 4999 simulated point patterns under 
the fitted Thomas process. 


in Section 3-3 we need a finite representation of which we treat as a latent process. This 
section considers a Bayesian Markov chain Monte Carlo missing data approach for the missing 
data tbs' and the unknown parameters pi > 0, // G k > 0, a > 0, and > 0 . 

Imagining that also a realization <1>5 = {(yi, rti),..., {yk,Uk)} had been observed, we detail 

below the calculation of the likelihood /[pL, P, tt, ct, cr^ | {xi,..., x^}, {(yi, ui),..., (y^, Ufc)}]. 

For the parameters, we assume independent prior densities p{pL),p{p),p{n),p{a),p{a‘^)', fur¬ 
ther prior specifications are given below. Hence the posterior density is 

p[pL,p, ft;,a,cj^{(yi,'ui),..., {yk,Uk)] \ {xi,.. .,Xn}] 

oc 1[PL, p, K, a, cr^ I {xi,..., x^}, {(yi, ui),..., (y^, Uk)}]p{pL)p{p)p{K)p{a)p{a‘^). (18) 

As a first ingredient of the likelihood, using the approximation <l>s' of <1>, we also approximate 
Xw by a finite Cox process Xw^s with driving random intensity function given by (10) with 
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k^±{-) = /(• I CT^). Conditional on $5, Xw^s is absolutely continuous with respect to the unit 
rate Poisson process on W, with density 


f[{xi,. .. ,Xn} I 4>s,a,fT^] =exp||lF| - j Kw{x \ 4>s, a, cr^) dxj | 4>s',a,o-^) 

(19) 

for finite point configurations {xi,... ,Xn} C W. 

For the second ingredient of the likelihood, notice that the distribution of $5 is absolutely 
continuous with respect to the distribution of a natural reference process <l>o,5 defined as fhe 
Poisson process on S wifh infensify funcfion = |u(‘^^|r((i/2)/(27r'^/^) wifh respecf fo 

fhe measure A(dy)r'rf_i(du), cf. Secfion 3-3. This reference process corresponds fo fhe case 
of an isofropic Poisson line process wifh unif intensify. The densify of <^5 wifh respecf fo fhe 
disfribufion of <l>o,s' is 


I PL,y,K\ 


= exp 


Us 


{xo{y,u) -x{y,u I PL, y, k)} X{dy)i^d-i{du) 


A x{yj,uj I pL,y,s) 
xo{yj,uj) 


for finite poinf configurafions {(yi, ui),..., {yk,Uk)} C S. Thai is, using fhe nofafion in Sec¬ 
fion 3-3, 


f[{{yi,ui),...,{yk,uk)} I PL,y,K\ 

A f 27r'^/2 

exp{-pL/(y,K)} n I e Ju 


(X 


where we have omiffed a consfanf nof depending on fhe paramefers. 
Combining (19)-(20) we obfain fhe approximafe likelihood 

1[PL,P, K, a, cr^ I {xi,..., Xn}, {{yi,ui), {yk,Uk)}] 


exp ||PF| - j Kw{x I 4>5,a,cj^) dxj JJ | 4's',a,cr^) 

X exY>{-pLl{y,s)} n I ^ Juj) j ■ 


( 20 ) 


( 21 ) 


Inserfing fhis info (18), we nofice fhaf fhe posferior densify is analyfically infracfable. A hybrid 
Markov chain Monfe Carlo algorifhm or Mefropolis wifhin Gibbs algorifhm, see e.g. Gilks el al. 
(1996), for posferior simulalions is proposed in Appendix C. Briefly, fhe algorifhm allernales be- 
Iween updaling each of fhe parameters and fhe line process, using a birfh-dealh-move Mefropolis 
Haslings algorifhm for fhe line process. 

To illusfrale fhe Bayesian approach we consider fhe Iwo-dimensional chapel dala sef in fhe lefl 
panel of Figure 2, using a uniform prior for bolh y = (cos (p, sin (p) and cr^, and Hal conjugated 
gamma priors for pi and a, see Figure 6. Our posferior resulls for pL, p, and a were sensilive fo 
fhe choice of prior disfribufion for n. For small values of n, i.e., values less lhan 30, meaningless 
posterior resulls appeared, since p was approximately uniform, and for p close fo zero, pi fended 
fo zero and hence a fended fo infinity. On fhe olher hand, very large values of k caused a very 
concenlrafed posferior disfribufion for p. As a compromise, after some experimenlalion, we fixed 
s = 40. 
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Figure 6. Prior and posterior distributions: The first four 
panels show the unnormalized prior density (solid line) and 
histogram for the posterior distribution of pL, ‘p, Q, and 
(T^, respectively. The final panel shows the histogram for 
the posterior distribution of p. 


For this model an extension of the observation window 14^ = [—.5, .5]^ to Wext = 
[—0.55,0.55]^ seemed large enough to aeeount for edge effeets. For the posterior simulations 
we used 200,000 iterations, where one iteration eonsists of updating all the parameters and the 
missing data. We eonsidered traee plots, whieh have been omitted here, for the parameters and 
information about the missing data, indieating that a burn-in of 5000 iterations is suffieient. 

There is a elear distinetion between the simulated posterior results for the parameters and 
the priors, ef. the first four panels in Figure 6. The posterior mean of ip (115.02°) is in elose 
agreement with the result of 117° found in Seetion 2-4, and a is unlikely to be larger than 0.02, 
indieating that the points are rather elose to the lines and that the ehoiee of IFext makes sense. 
The final panel in Figure 6 shows a good agreemenf befween fhe number of ehapels (110) and 
fhe posferior mean of fhe infensify p (103.7), fhough fhere is some uneerfainfy in fhe posferior 
disfribufion of p. Moreover, fhe posferior means of and a are 12.9 and 8.4, respeefively, whieh 
eombined wifh (15) resulf in fhe esfimafe 108.4 for p. 

To illusfrafe fhe usefulness of fhe Bayesian mefhod in defeefing linear sfruefures. Figure 7 
shows a posterior kernel esfimafe of fhe densify of lines wifhin W. The esfimafe visualizes where 
fhe hidden lines eould be, i.e., fhe lighfer areas, and overall fhey agree wifh fhe poinf paffern of 
ehapels, whieh is superimposed in fhe figure, fhough in fhe upper righf eorner of fhe observation 
window fhere is some doubf abouf whefher fhere should be a single or fwo elusfers of poinfs. 
Speeifieally, fhe esfimafe is obfained from 100 posferior iferafions wifh an equal spaeing, and if 
is fhe average of binary pixel represenfafions of fhe line proeess, where a pixel has value 1 if if is 
inferseefed by a line, and value 0 ofherwise. 


4. Discussion 

4-1. Choice of structuring element 

Ripley’s TT-funefion has a ball as sfruefuring elemenf buf fhis is nof useful for defeefing 
anisofropy. Direefional iT-funefions have been suggested using a seefor annulus (Ohser & 
Sfoyan, 1981) or a double eone (Redenbaeh ef ah, 2009) as fhe sfruefuring elemenf, while we 
have suggesfed a eylinder. Anofher suggesfion eould be an ellipsoid. 

A defailed eomparison of fhe eylindrieal AT-funefion wifh fhe iT-funefion in Redenbaeh ef al. 
(2009) using a double eone as fhe sfruefuring elemenf is given in a feehnieal reporf by F. Safav- 
imanesh and C. Redenbaeh from 2016. They eonelude fhaf in sifuafions where fhe anisofropy 
is pronouneed, a good ehoiee of sfruefuring elemenf may be imporfanf and will depend on fhe 
applieafion af hand. In ease of geomefrie anisofropy (Mpller & Toffaker, 2014), an ellipsoid is 
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Figure 7. Posterior kernel estimate of the density of lines. 
For comparison, the chapel point pattern data set is super¬ 
imposed. 


appropriate; when there is a eolumnar strueture as under the Poisson line eluster point proeess 
model or in the data sets eonsidered in this paper, an elongated eylinder is appropriate; while if 
regular point proeess models are eompressed, then a double eone is appropriate. F. Safavimanesh 
and C. Redenbaeh emphasize the importanee of an appropriate ehoiee of the seale and shape pa¬ 
rameters used to speeify the strueturing element, i.e. r and t in ease of Ku{r,t). They notiee 
that prior information, e.g., the diameter of the clusters/mini-columns of points in ease of the 
pyramidal eells (Rafati et ah, 2016), ean be used to determine interesting ranges of r values. 


4-2. Non-stationary case 

In some applieations it is relevant to use a non-eonstant intensity funetion p{x). Suppose 
we assume seeond order intensity reweighted stationarity (Baddeley, Mpller & Waagepetersen, 
2000) and g{x) still denotes the pair eorrelation funetion. This means intuitively, that if xi,X 2 G 
IR^ are distinet loeations and Bi, B 2 are infinitesimally small sets of volumes dxi, dx 2 and eon- 
taining xi, X 2 , respeetively, then p{xi)p{x 2 )g{xi — X 2 ) dxi dx 2 is the probability for X having 
a point in eaeh of Bi and i? 2 - Our definition (1) still applies, while (2) beeomes 


Ku{r,t) 



E 

xi,X2eX: xi^a;2 


l{xi £ PF, X 2 - Xi £ Cujr, t)} 

p{xi)p{x2) 


r > 0, t > 0, 


whieh in turn ean be used when deriving non-parametrie estimates. 

Assumption (B2) may be relaxed to obtain a non-stationary Poisson line eluster point proeess 
model for X, assuming that for eaeh line the Poisson proeess Qi has intensity funetion ai(x) = 
a(x) for X G /*, where a is a non-negative funetion whieh is loeally integrable on any line in IR*^. 
Then (6) should be replaeed by 


00 

A(x) = y^a{(x ■ Ui)ui + gi}k^±{p^±(x - Pi)}, x G (22) 

^ ^ I i 

i=l 


However, this non-stationary extension of the model will be harder to analyze, e.g., moment 
results as established in Seetion 3-4 for the stationary ease will in general not easily extend. 
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and it turns out that the model is not seeond order intensity reweighted stationary exeept in a 
speeial ease diseussed below. Moreover, while our statistieal methodology in Seetion 34 ean be 
straightforwardly extended, the Bayesian eomputations in Seetion 3-6 beeome harder. 

In (22) assume that the intensity funetion a{y^^\ ... does not de¬ 
pend on y = ..., 0) G H. Let c = fj a{x^'^^} dx^^\ Then, using a notation as in 

Seetion 34, it ean be shown that X is seeond order intensity reweighted stationary, Xj is still a 
Neyman-Seott proeess, while the x^^'^^-eoordinates of the points in Xw are independent and iden- 
tieally distributed with density and they are independent of Xw- Therefore statistieal 

inferenee simply splits into modelling the density a{x^'^^}/c based on the x^^'l -eoordinates of the 
points in Xw and inferring {c,pL,o'^) by eonsidering Xj along similar lines as in Seetion 34 
but with replaeed by c. 
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Appendix A: Intensity and rose of direction for a Poisson line process 


First we give the definition of the intensity pi and the rose of direetion TZ for a general sta¬ 
tionary line proeess L = {( 1 ,( 2 , • •.} in IR*^. Let | • |i denote one-dimensional Lebesgue measure, 
dt Lebesgue measure on the real line, 4 C IR'^ an arbitrary Borel set with volume |4| G (0, 00 ), 
and an B c S"*-! arbitrary Borel set. Then by definition and sinee L is stationarity. 


PL = E 


' 00 

E 

\i=l 


|fin4|i/|4| 


(23) 


does not depend on the ehoiee of A, and provided 0 < < 00 , 


77(B) = E 


\li n 4|il(tti G B)/ (pl|4|) 


(24) 


does not depend on the ehoiee of A and is seen to be a probability measure. 

Seeond we assume that L is a stationary Poisson line proeess as in Seetion 3-1. Then 


B|f;|f,n4|il(nG B)J 


J l{yi + tUi G A, Ui G B) dtj 

(25) 


^4 

^ J J ^{y + tu G A, u G B) dt X{dy) M(du) 

(26) 


= /3\A\ [ l/|n(‘^)|M(du). 

Jb 

(27) 


Here (25) follows from the phase representation of L (see Seetion 34), (26) from the Slivnyak- 
Meeke theorem for the Poisson proeess (see e.g. Mpller & Waagepetersen (2004)), and (27) 
sinee is the Jaeobian of the mapping (t, y) y + tu with (t, y) G IR x 77. When B = 

we obtain from (23) and (27) the first equation in (7). This together with 0 < /3 < 00 and 
0 < J^d-i M(dtt) < 00 imply that 0 < pL < 00 . Thereby the seeond first equation in 

(7) follows for any Borel set B C 
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Appendix B: Moment results for a Poisson line cluster point process 
For Section 34 it remains to verify (15)-(16). 

Proof of (15); By (6), (14), and the Slivnyak-Mecke theorem for the Poisson process, 

P = a/3j j ku{-Puiy)} K'iy) M{du). (28) 

Let Id be the d x d identity matrix, and Od-i the origin in For u G let A{u) = vv^ 
where v is the subvector consisting of the first d — \ coordinates of u and ^ denotes transpose of 
a vector or matrix. The Jacobian of the linear transformation 

Pu{y) = {Id - uu^) y, y ^ H, 

is the square root of the determinant of the (d—1) x (d—1) matrix 

Q = [{Id- {Id-1 Od-l)^} {Id - UU^) {Id-1 Od-l)^ 

= {Id-1 Od-l) {id — UU^) {Id-1 Od-lf" 

= Id -1 - A{u). 

Since A{u) is symmetric of rank at most one and has trace tr{4(ri)} = + ... + 

{.y(d-i )}2 _ determinant of Q is 1 — tr{4(u)} = Combining this with 

(28) we obtain 

p = a/3j j ku{-y)/\u^‘^^\Xu{dy) M{du). 

Thereby (15) easily follows from the first identity in (7). 

Proof of (16); By (6) and (14), 


p^g{x) = a^E 


'^K^{Pu^{-yi)}K^{Pu^{x - Vj)} 

< ^ XX 3 3 

i¥=j 


+ a^E 


'^K^{Pu^{-yi)}K^{Pu^{x - yi)} 

< ^ XX XX 


= p^ + a^l3 j k^±{p^±{-y)}k^^{p^±{x - y)} X{dy) M{du) 


(29) 


using the extended Slivnyak-Mecke theorem for the Poisson process and the proof of (15) to 
obtain that the first expectation is equal to p^, and the Slivnyak-Mecke theorem for the Poisson 
process to obtain that the second expectation is equal to the last term. Combining (15) and (29) 
with the result for the Jacobian considered above, we obtain (16). 


Appendix C: Hybrid Markov chain Monte Carlo algorithm 
This appendix details the Markov chain Monte Carlo algorithm for the Bayesian approach 
considered in Section 3-6, with independent prior densities for the parameters and posterior den¬ 
sity given by (18)-(21). As in Section 3-6, we consider conjugated gamma densities p{a) and 
p{pl), and denote their shape parameters by oi and 02 and their inverse scale parameters by bi 
and 62 , respectively. The remaining parameters have no (well-known) conjugate priors, cf. (19) 
and ( 20 ), and thus we consider generic prior densities p{p.), p{k), and p{o'^). 
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In each iteration of the Markov chain Monte Carlo algorithm we update first each of the pa¬ 
rameters and second the missing data. We use a Gibbs update for a respective p^, noting that 
the conditional distribution of a given the rest is a gamma distribution with shape parameter 
ai -I- n and inverse scale parameter bi -\- fiPu-!- ~ hj) I and the conditional 

distribution of pL given the rest is a gamma distribution with shape parameter 02 + k and inverse 
scale parameter 62 + -^ (p, k)- Below we describe the individual proposals and Hastings ratios for 
the remaining parameters and the missing data (in the case of Section 3-6 where the value of k 
is fixed, we can of course jusf ignore fhe updafe of k described below). As usual, for each fype 
of updafe, fhe proposal is accepfed wifh probabilify min{l,i?}, where R is fhe corresponding 
Hastings ratio. We denofe {pL, p, k, a, ..., {yk,Uk)}) the current state of the al¬ 

gorithm, where n > 1, k > 1 (since A: = 0 implies n = 0, which is not a case of interest), and 

tii) n Vkext 7^ 0, f = 1, •. •, (c. 

For p, K, and a‘^, we use Metropolis random walk updates, with a von Mises-Fisher pro¬ 
posal p' ~ /(• I p,Ko) and normal proposals k' ~ A^(k, cjg^) and ~ A^(it^,cJq 2 ), where 
Ko, Cq 1 ) ^0 2 > 0 are tuned so that the mean acceptance probabilities are between 20-45% (as 
recommended in Roberts, Gelman & Gilks (1997)). The Hastings ratios for the acceptance prob¬ 
abilities are 

k 

D P{P') r . j, , ^llTjf{Uj\p',K) 


Rk - !(«' > exp 


f{Uj I p, k') 


[pL {liP, -) - HP, -')}] n I ’ 


and 


R ^2 = i(a ^ > 0) ^^^2 exp ( a ^ [ f{p^±ix - yj) \ aHdx 

FW ) y ^ Uw ^ 




- / f{Pu^ix-yj) \(J }dx 

Jw ^ 


n 


EUf{PuHx^-yH\<^'H 


EU HPuHx^ - yj) \ 


For R ^2 each integral is calculated by a simple Monte Carlo method after making a change of 
variables from x to its Cartesian coordinates in a system centered at yj and with axes given by 
Uj and uj-. 

For the missing data, we adapt the birth-death-move Metropolis-Hastings algorithm in Geyer 
& Mpller (1994) as follows. Each of the birth/death/move proposals happens with probability 
1/3 and consists of the following action. A birth proposal is the proposal of adding a new point 
{y, u), where u ~ /(• | p, k) and y conditional on u is uniformly distributed on J^. Then, as 
explained below, the Hastings ratio is 


i?birth = u) n Wext / 0} 


X exp 

-« / f{Pux{x - y) \ crHdx 

n 

n 

f{Puxixi-y)\(jH 


Jw 

i=l 

Ei=i f{Pux{xi-yj) 1 (xH 

. J 


(30) 
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To stress the dependenee on {{yi,ui),... ,{yk,Uk)] and {y,u), write i?birth = 
^birth(yi, • • •, yfc, Uk] y, u) (obviously, it also depends on pL, a, and {xi,..., Xn}). 
A death proposal is the proposal of generating a uniform j G {1,..., A:} (provided A: > 1; if 
A = 1, we do nothing and keep the eurrent state) and deleting {yj,Uj). Then, as explained below, 
the Hastings ratio is 

A^death — 1/A?birth (?/l) ^1) • • • ; 2/j —1 > —1) Vj+l^ • ’ ' ) t Vj ; Uj ). (31) 

Finally, a move proposal is the proposal of seleeting a uniform j G {1,..., A} and replaeing 
{yj,Uj) by {y'pu'j), where n' ~ /(• | ;U,k) and y'- eonditional on u' is uniformly distributed 
on J^i. Sinee this ean be eonsidered as first a death proposal and seeond a birth proposal, the 
Hastings ratio is 

^ A?birth(yi) 'U'l) • • • ) Vj—li 1) 2/i+l) ^i+1 • • • ) Vki Uk] Vj, Uj) 

-^move — 7 T • 

ribirth(2/l; 'W-l > • • • jUj—li'^j—liUj+li '^j+1 ■ ■ • jVki Uj > '^j ) 

It remains to explain how we obtained the Hastings ratios (30)-(31), where we notiee the 
following faets. The referenee Poisson proeess ‘ho,s’ has intensity measure 

C{dy,du) = |u("*)|^^^A(dy)jVrf_i(du). 

Further, eonditional on the data Xw = {xi,... ,Xn} and the parameters p^, p, k, a, the 
target proeess <I> 5 ' has density 

f[{{yi,Ul), ..., {yk,Uk)} I {xi,.. .,Xn},a,a‘^,pL,p,K] 

oc /[{xi,... ,Xn} I {{yi,ui),... ,{yk,Uk)},a,a‘^]f[{{yi,ui),... ,{yk,Uk)} \ Pl,^,^] 

with respeet to the distribution of <l>o, 5 . Furthermore, if a birth {y, u) is proposed, then it has 
density 

fju I fi,K)Hy€Ju)/XiJu) 

’ \uW\r{d/2)/{2Tr<^/^) 

with respeet to C- Consequently, by Geyer & Mpller (1994), the Hastings ratio for the proposed 
birth is 

^ f[{{yi,ui), ■ ■ ■, {yk,Uk), {y,u)} I {xi, ■ ■. ,Xn},a,cr‘^,PL,^k,K] ^ 1/(A + 1) 

f[{{yi,Ul),---,{yk,Uk)}\{xi,...,Xn},a,a‘^,PL,P,K] f{y,u\p,K) 


2'k^I^ 1 

T(d/2) {k + l)f{y,u\ p,k) 


X exp 

-« / f{Pu^{x-y)\a^}dx 

n 

n 

f{Pux{xi-y)\(T‘^} 


Jw 

i=l 

Ej=i f{Pufi^i-yj) 1 


whieh is equal to (30). Thereby, refering again to Geyer & Mpller (1994), we obtain (31). 
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